        subroutine UHF(NORB,F_a,F_b,P_a,P_b,n_el_a,n_el_b,Num_ele,S,over_m,Hcore,X_T,U,E_nu)
        implicit none
        integer :: i,j,k,l,count_iter
        integer :: NORB,Num_ele,n_el_a,n_el_b
        integer :: n_occ_a(NORB,NORB),n_occ_b(NORB,NORB)
        real*8  :: X_T(NORB,NORB),X_A(NORB,NORB)
        real*8  :: Hcore(NORB,NORB),over_m(NORB,NORB)
        real*8  :: U(NORB,NORB,NORB,NORB)  !!!two-ele inte
        real*8  :: F_a(NORB,NORB),F_b(NORB,NORB)
        real*8  :: F_aa(NORB,NORB),F_bb(NORB,NORB)
        real*8  :: P_a(NORB,NORB),P_b(NORB,NORB)
        real*8  :: P_a1(NORB,NORB),P_b1(NORB,NORB)
        real*8  :: C_a(NORB,NORB),C_b(NORB,NORB)
        real*8  :: eig_a(NORB,NORB),eig_b(NORB,NORB)
        real*8  :: tem(NORB,NORB),temp(NORB,NORB)
        real*8  :: orb_ae(NORB),orb_be(NORB),Sem12(NORB,NORB)
        real*8  :: Ptolold(NORB,NORB),Ptolnew(NORB,NORB)
        real*8  :: E_ele1,E_ele,delta,delta1,E_nu,rms
        logical :: logical_e 
!!!===
!!!! spin contamination
        real*8  :: S_2_exact,S_2_UHF
        integer :: S  !!!! spin multi
        
        call set_nel(NORB,S,Num_ele,n_el_a,n_el_b)
        write(21,*)'Alpha ele',n_el_a   !n_el_a=9;
        write(21,*)'Beta ele ',n_el_b   !n_el_b=7;
!!!inital guess
        write(21,*)'****UHF program activate****'
        !Form density martix
        !Use Hcore as initial guess for Hamiltonian
        !********Alpha*****************************
        X_A=transpose(X_T)    !!!X_T=U*S^(1/2)
        F_a=Hcore
        F_aa=matmul(X_A,matmul(F_a,X_T))
        tem=0.0d0
        tem=F_aa
        !write(*,*)tem
        call DIAG(tem,eig_a,NORB)
        C_a=matmul(X_T,tem)
        call Density(C_a,P_a,NORB,n_el_a)
        !**********end Alpha**********************
        !**********Beta*************************** 
        do i=1,NORB
                do j=1,i
                if(i==j)then
                F_b(i,j)=Hcore(i,j)
                else
                F_b(i,j)=0.5d0*1.7d05*over_m(i,j)*(Hcore(i,i)+Hcore(j,j))
                F_b(j,i)=F_b(i,j)
                endif
                enddo
        enddo
        F_bb=matmul(X_A,matmul(F_b,X_T))
        tem=0.0d0
        tem=F_bb
        call DIAG(tem,eig_b,NORB)
        C_b=matmul(X_T,tem)
        call Density(C_b,P_b,NORB,n_el_b)
        call E_UHF(NORB,F_a,F_b,P_a,P_b,Hcore,E_ele1)
        !**********end Beta**********************
        !!!end inital guess
        logical_e=.true.
!!!cyc
        count_iter=0
        do while(logical_e)
        count_iter=count_iter+1
        F_a=Hcore;F_b=Hcore
         do i=1,NORB
           do j=1,NORB
              do k=1,NORB
                 do l=1,NORB
                    F_a(i,j)=F_a(i,j)+(P_a(k,l)+P_b(k,l))*U(i,j,k,l)-P_a(k,l)*U(i,k,j,l)
                    F_b(i,j)=F_b(i,j)+(P_a(k,l)+P_b(k,l))*U(i,j,k,l)-P_b(k,l)*U(i,k,j,l)
                 enddo
              enddo
           enddo
         enddo
        !!!**********Alpha********************
        F_aa=matmul(X_A,matmul(F_a,X_T))
        tem=0.0d0
        tem=F_aa
        call DIAG(tem,eig_a,NORB)
        C_a=matmul(X_T,tem)
        call Density(C_a,P_a1,NORB,n_el_a)
        !!!**********end**********************
        !!!**********Beta*********************
        F_bb=matmul(X_A,matmul(F_b,X_T))
        tem=0.0d0
        tem=F_bb
        call DIAG(tem,eig_b,NORB)
        C_b=matmul(X_T,tem)
        call Density(C_b,P_b1,NORB,n_el_b)
        !!!**********end Beta*****************
        call E_UHF(NORB,F_a,F_b,P_a1,P_b1,Hcore,E_ele)
        call delta_P(NORB,P_a,P_a1,delta)
        call delta_P(NORB,P_b,P_b1,delta1)
        !call delta_P(NORB,p_a+p_b,p_a1+p_b1,delta) 
        ptolold=P_a+P_b
        ptolnew=P_a1+P_b1
        rms=0.0d0
        do i=1,NORB
                do j=1,NORB
                rms=rms+sqrt((Ptolnew(i,j)-Ptolold(i,j))**2)
                enddo
        enddo
        rms=rms/NORB
        !if(abs(E_ele-E_ele1)<1.0e-6.and.(delta<1.0e-6).and.(delta1<1.0e-6))then
        if(abs(E_ele-E_ele1)<1.0e-8.and.rms<1.0e-5)then
        logical_e=.FALSE.
        write(21,*)'****************************'
        write(21,*)'UHF cycle converged!'
        write(21,*)'Number of iterations:',count_iter
        write(21,*)'Electronic energy   :',E_ele
        write(21,*)'Nuclear energy      :',E_nu
        write(21,*)'UHF Energy          :',E_ele+E_nu
        write(21,*)'***************************'
        write(21,*)
        write(21,*)"Orbital Energy:      "
       !call DIAG_1(F_aa,orb_ae,NORB)
       !call DIAG_1(F_bb,orb_be,NORB) 
       !do i=1,NORB
       !write(21,'(2X,A3,I2,A1,3X,2F16.12)')"orb",i,":",orb_ae(i),orb_be(i)
       !enddo
       !write(*,*)'1111'
        endif 
        E_ele1=E_ele
        P_a=P_a1;P_b=P_b1
        
        enddo

        !*********************************************************
        !******************* SPIN CONTAMINATION ******************
        !*********************************************************

        !> Calculate the expectation value of S^2 
        S_2_exact=(n_el_a - n_el_b)/2.0d0 *&
                           ((n_el_a - n_el_b )/2.0d0 + 1.0d0)
        S_2_UHF=0.0d0
        temp=matmul(matmul(transpose(C_a),over_m),C_b)
        do i=1,n_el_a
           do j=1,n_el_b
                    S_2_UHF=S_2_UHF+temp(i,j)**2.0d0
           enddo
        enddo
        S_2_UHF=S_2_exact+n_el_b-S_2_UHF
        write(21, "(A)") "Spin contamination:"
        write(21, '(A, F20.14)') '  <S>^2                   :',S_2_UHF
        write(21, '(A, F20.14)') '  <S>^2 expectation value :',S_2_exact
        write(21, '(A, F20.14)') '  Spin contamination      :',abs(S_2_UHF-S_2_exact)
        write(21,*)
        endsubroutine UHF
         
 
        subroutine Density(C_new,P_new,NORB,Num)
        !input
        implicit none
        integer :: NORB,Num
        integer :: i,j,k,x
        real*8  :: C_new(NORB,NORB)      !New coefficient matrix
        !output
        real*8 :: P_new(NORB,NORB)       !New density matrix
        P_new=0.0d0
        do i=1,NORB
           do j=1,NORB
              do k=1,Num
                 P_new(i,j)=P_new(i,j)+C_new(i,k)*C_new(j,k)
              enddo
           enddo
        enddo
        endsubroutine Density

        subroutine E_UHF(NORB,F_a,F_b,P_a,P_b,Hcore,E_ele)
        integer :: i,j,k,l,NORB
        real*8  :: F_a(NORB,NORB),F_b(NORB,NORB)
        real*8  :: P_a(NORB,NORB),P_b(NORB,NORB)
        real*8  :: Hcore(NORB,NORB),P_t(NORB,NORB)
        real*8  :: E_ele
        P_t=0.0d0
        P_t=P_a+P_b
        E_ele=0.0d0
        do i=1,NORB
           do j=1,NORB
              E_ele=E_ele+0.5d0*(P_t(i,j)*Hcore(i,j)+P_a(i,j)*F_a(i,j)+P_b(i,j)*F_b(i,j))
           enddo
        enddo
        
        endsubroutine E_UHF 

        subroutine set_nel(NORB,S,Num_ele,n_el_a,n_el_b)
        implicit none
        integer :: NORB,S !S:spin multiplicity
        integer :: Num_ele,n_el_a,n_el_b
        if(S==3.and.mod(Num_ele,2)==0)then
                n_el_a=Num_ele/2+1
                n_el_b=Num_ele/2-1
        else if(S==2.and.mod(Num_ele,2)==1)then
                n_el_a=(Num_ele-1)/2+1
                n_el_b=(Num_ele-1)/2
        endif

        endsubroutine set_nel

